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Abstract. Wc study local moment formation for adatoms on bilayer graphene 
(BLG) within a mean field theory of the Anderson impurity model. The 
wavefunctions of the BLG electrons induce strong particle-hole asymmetry and 
band dependence of the hybridization, which is shown to result in unusual features 
in the impurity model phase diagram. We also study the effect of varying 
the chemical potential, as well as varying an electric field perpendicular to the 
bilayer; the latter modifies the density of states of electrons in BLG, and, more 
significantly, shifts the impurity energy. We show that this leads to regimes in the 
impurity phase diagram where local moments can be turned on or off by applying 
modest external electric fields. Finally, we show that the RKKY interaction 
between local moments can be varied by tuning of the chemical potential (as 
has also been suggested in monolayer graphene) or, more interestingly, by tuning 
the electric field so that it induces changes in the band structure of BLG. 
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1. Introduction 

Single layer graphene hosts a plethora of phenomena that arise from the Dirac- 
like band dispersion and chirality of its low-energy quasiparticle excitations [U 
It is interesting to explore how these unusual single particle properties impact the 
physics of adatoms on graphene. The combination of adatom-graphene hybridization 
and Hubbard-like interactions on the adatom has been studied in the context of 
local moment formation [3j H], Kondo physics El [7j [9], RKKY interactions 
HDim D21H31 [HH31HB], and adatom Positional ordering [TO OH QUI [SO] . The study 
of adatoms on monolayer graphene is also of interest to the nanoscience and quantum 
computation communities given the possibility to control local moment physics, and 
adatom-adatom spin and density interactions, by varying the carrier concentration via 
gating [2"T] . 

In contrast to monolayer graphene, bilayer graphene (BLG), which has Bernal 
stacking of single layers, has an extra tuning parameter. Using a dual-gate geometry, 
shown in figure[T^, enables one to separately tune the chemical potential and an electric 
field perpendicular to the layers, which is equivalent to separately tuning the potential 
on each of the two layers of BLG. While tuning the chemical potential modifies the 
carrier concentration, applying an electric field normal to the layers generates a gap 
in the band structure of BLG [21 OH [Ml \M (23 (Figure [J). (We will refer 
to the potential difference between the two layers, induced by this electric field, as 
the 'bias'.) Such a tunable gap system enables one to envision device applications 
and the ability to dynamically control various states in BLG [HI HH] This tunability 
also allows for the study of interesting fundamental physics — for instance, it has 
been shown that engineering the electric field to flip direction (from pointing up to 
pointing down) as a function of position leads to localized one-dimensional modes at 
the kink in the bias [301 Ell GUI ■ We have shown in recent work that incorporating 
interaction effects converts this 'nanowire' into a 2-band Tomonaga-Luttinger liquid 
whose properties, such as Luttinger parameters and mode velocities, can be controlled 
by the bias strength [33] . 

In this paper, we study adatoms in BLG. We examine local moment formation on 
the adatoms, RKKY interaction between such local moments, and how these effects 
can be controlled by tuning the chemical potential and a applying perpendicular 
electric field. 

Our work goes beyond Ref . |34] . which studied local moment formation for site- 
centered adatoms on BLG, in several important respects, (i) We consider adatoms that 
are positioned at the center of a hexagonal plaquette on one of the layers. The study 
of this configuration is motivated by a recent ab initio study of adatoms in monolayer 
graphene that indicates plaquette centered impurities are generally more energetically 
favourable than on-site impurities [35] . We expect a similar situation to hold in BLG. 
(ii) An applied electric field is shown to directly tune the impurity energy. This is 
because an impurity position will, in general, be located closer to the top layer of BLG. 
Accounting for this impurity energy shift allows us to identify regions of the phase 
diagram where local moment formation can be turned on and/or off by the application 
of a perpendicular electric field, (iii) For a particular impurity level chosen so that its 
renormalized (with self-energy corrections) energy level lies in the middle of gap in 
presence of the bias, we construct phase diagrams at zero, positive and negative bias 
by sweeping the chemical potential. The resulting phase diagram exhibits the onset 
of a Coulomb-blockade phase where any arbitrarily small U results in the formation 
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Figure 1. (a) Bilayer graphene in a dual-gate configuration, (b) Schematic 
diagram of a plaquette-centered (large, red) adatom impurity on the top layer of 
bilayer graphene. (c) Cross-section of the dispersion relation for unbiased (Left) 
and biased (Right) graphene close along k y = through the K-point (A = and 
A = 0.025t, respectively). Inset: The two unique K-points and the cross-sectional 
cut are indicated in the Brillioun zone. 

of local moments, (iv) As a consequence of the chiral wavefunctions of BLG and the 
fact that the plaquette centered impurity adatom couples to many sites, the coupling 
between the impurity and the quasiparticles of BLG has strong momentum and band 
dependence. This affects many of the details of the phase diagram. For instance, the 
self-energy develops a large real part that has nontrivial frequency dependence, and 
substantially renormalizes the position of the impurity spectral peak in a manner that 
depends on the chemical potential and the applied bias. We provide detailed a physical 
explanation for how this affects the resulting phase diagrams, which were not provided 
in reference [33]. Furthermore, to better illustrate the effect of the wavefunctions and 
chirality of BLG on the phase diagrams, the BLG system is contrasted with a fictitious 
system of non-chiral fermions with the same DOS and dispersion relation, (v) We go 
beyond the issue of local moment formation to address the tunable RKKY interactions 
between such local moments on BLG. 

We begin, in Section II, by introducing the Anderson impurity model specific to 
BLG. Section III summarizes the Anderson mean field theory formalism. Armed with 
this background, in Section IV we construct the impurity model phase diagrams for 
plaquette-centered adatoms (shown schematically in figure [TJj) . To highlight some 
of the unusual features of these phase diagrams, we contrast it with an impurity 
model of a fictitious system of electrons that have an identical dispersion but a band- 
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independent coupling to the adatom. Finally, in Section V, we discuss the RKKY 
interaction, and its tunability, for local moments on BLG. 

2. Adatom model in bilayer graphene 

Consider an adatom on BLG, described by the Anderson impurity model [36], 

#BLG = ( £ ks - M)cLa C k S a> (!) 

k,S,fT 

H imp = ^2{e d - n)dld a + Un dt n dV (2) 

(X 

H mix = -^Xr( c L^+4 c r<r)- (3) 

Here e ks is the BLG electron dispersion for electrons labelled by momentum k and 
band index s. We assume a minimal model for the BLG dispersion that includes a 
nearest-neighbor hopping amplitude, t, to sites on the same layer, and an intcrlaycr 
hopping amplitude, t±, between the two sites that sit one on top of the other. 
Henceforth, we set t — 1 and note that t ps 3 eV and t±/t s» 0.15 in BLG. In H[ nip , 
we denote the impurity energy by e^, while U denotes the electron-electron repulsion 
on the impurity site. BLG electrons at sites r can hop on or off the adatom impurity 
with an amplitude \ r - We assume a common equilibrium chemical potential fx for the 
impurity and BLG electrons. The complete Hamiltonian for unbiased BLG is then 
given by H = H B lg + H imp + H mix . 

Electronic structure studies of transition metal adatoms on monolayer graphene 
suggest that the low-energy configuration of many types of impurities corresponds to 
the adatom residing at the center of a hexagonal plaquette [35] . We therefore fix the 
adatom position to be at the plaquette center on the top layer (labelled £ — 1) of 
BLG, as shown in the schematic diagram on the right in figure [l] For simplicity, we 
assume that Xr~X f° r the set of sites {r„}, which includes the six nearest neighbor 
plaquette sites in layer- 1 and the site on layer-2 that lies directly below the adatom, 
and Xr = for all other sites. This simplifying assumption about the impurity model 
allows us to focus on unconventional features of local moment formation intrinsic to 
bilayer graphene. Future density functional studies would be useful in incorporating 
details of the impurity atomic orbitals. Turning to the mixing Hamiltonian H m [ K 
which allows the impurity electrons to hybridize with the BLG electrons, let us set 

V ks ^X £ ks (r), (4) 

r={r„} 

where </> ks (r) denotes the wave function at site r for electrons in band-s and momentum 
k. We then obtain 

k.s,cr 

+ Y.( e d~ri d U* + Un dt n di . (5) 

a 

While the impurity model Hamiltonian in BLG looks similar to that in conventional 
systems or monolayer graphene, there are two important new ingredients in the 
impurity physics of BLG with plaquette centered impurities. 
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Figure 2. Coupling of the impurity to the four bands (ordered from lowest 
to highest energy and scaled by system size), (a) |Vki|> (b) IVfol, (c) |Vk3|, and 
(d) |Vk4|, with impurity hopping strength \ = 0.3t. Dotted line indicates the 
Brillouin zone. 

First, for BLG (or multilayer graphene), as opposed to monolayer graphene, one 
can tune the density of states by applying an electric field perpendicular to the layers. 
Let A denote the potential imbalance between the top and bottom layer induced by 
the electric field. Assuming that the adatom is at the same height as the top layer, 
this leads to an extra term in the Anderson Hamiltonian 

ffbias =-f £ {-\fc\ ta C ria +^£44, (6) 

where denotes the sites in the top (i = 1) and bottom (£ — 0) layers. In writing this 
modification to the Hamiltonian, we have assumed that \ an d t remain unchanged 
in the presence of an electric field. If intercalation of the impurity occurs, this will 
reduce the shift in the impurity energy, but will always be nonzero on grounds of 
the crystal symmetry. Incorporating the bias in this way thus has three effects: (i) 
a renormalization of the BLG dispersion; (ii) a modification of the hybridization 
Vk s through a change in the BLG quasiparticle wavefunctions; and (iii) a shift the 
impurity energy to e d + A/2. We will refer to the renormalized BLG dispersion and 
the hybridization as e ks (A) and Vk s (A) respectively. It is well-known that such a 
bias in BLG can open a band gap and significantly change the low-energy density of 
states; what is perhaps not appreciated is that this also effectively tunes the impurity 
energy in multilayer graphene. The last term in equation [6] describing this effect was 
not present in reference [33] and it will be shown to have a remarkable effect on local 
moment formation in presence of a bias. 

A second important difference arises from the tunneling matrix elements, Vk s , for 
the four bands of the bilayer. As shown in figure [2] these matrix elements display 
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strong band- and momentum-dependence, which does not appear for the site-centered 
impurities discussed in reference |34) . The rich structure of the coupling between 
the chiral BLG quasiparticles and the impurity site leads to a number of differences 
in the impurity model phase diagram when compared with conventional non-chiral 
fermions with a similar density of states, where we simply replace </>ks(r) ~ exp(ik • r) 
in equation [4] 

3. Mean field theory 

A mean field treatment of the adatom impurity model is obtained, following Anderson 
[36] . by setting 

Un dt n di = U ^(\Pd- am d) n da ( 7 ) 

where p d = J2a( n dcr)^ an ^ m d — \ J2a a ( n da)- Let us tlien define 



Zd« =e d -fi + U(^-am d ) (8) 



rPd 
2 

4 5 (A) = e ks (A)- M . (9) 
With this mean field approximation, the entire Hamiltonian splits into two single 
particle impurity Hamiltonians, one for each spin, with 

#imp =(Zd«+^)dUa (10) 

#BLG= ££ ks (A)cL<W (11) 
k,s 

= -£(^(A)4 SCT d CT +y k * s (A)4c w ). (12) 

k 

These are coupled together by the self-consistency conditions that fix £ dcr via m d and 
p d . The single particle Green function for the impurity is given by 

G »n) = ~ ~ \ , (13) 

where the impurity self-energy is given by 

:(A)P 

We can analytically continue this to the real frequency axis by setting iw n — > lj + i0+ 
to obtain the real and imaginary parts of the self-energy T, d (uj). We can then compute 
at T = 

1 f° 

p d = / cL>£ Im G dd (iuj n ^uj + i0 + ), (15) 

^ J —00 „ 



x ^ |V- ks (A)| 2 

S rf ^») = ^ (14) 



1 T rf^.lm G5 d ( 



iw n ^w + i0 + ). (16) 

Within this mean field approach, the presence of a local moment on the impurity is 
signalled by a self-consistent solution with a nonzero m d . 

Alternatively, it is possible to self-consistently solve the mean field Hamiltonian 
using exact diagonalization for small system sizes. All of the phase diagrams in the 
next section were checked for consistency using this method. 
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4. Local moment formation 

Using the above mean field theory enables us to study local moment formation on an 
impurity atom residing on BLG. Since the BLG band structure can be tuned by the 
electric field, we choose to define Tq = Trx 2 /t as a rough scale for the impurity level 
broadening in the absence of interactions. Thus, Tq remains fixed for a given x even 
as the electric field and chemical potential are varied. In this section, we begin by 
discussing the case when A = (i.e. without an applied electric field perpendicular to 
the layers). Phase diagrams are constructed by varying t d and U for fixed \ = 0.3i 
(which implies \ ~ 1 e V in conventional units) with various choices of the chemical 
potential. Next, we consider how varying A can be used to tune the phase diagrams. 
After which, we discuss an alternative phase diagram for an impurity with a fixed 
bare energy level (although the actual energy level will be modified in the presence of 
a bias) with various choices of A. To construct these phase diagrams, \x and U are 
varied, while and \ (— 0.3i) are kept fixed. 

We have checked that varying \ modestly makes no qualitative changes to various 
features in the phase diagram, although it does shift the phase boundaries as expected. 
We ascribe the complexities of the impurity model phase diagram in BLG to the 
effective momentum- and band-dependent mixing Vk s . As we discuss below, the strong 
variation of this coupling between different bands results in particle-hole asymmetry 
of the impurity model phase diagram via the impurity self-energy. This is despite the 
fact that in the simplest tight-binding parameterization, which we have considered, 
the BLG band dispersion itself is particle-hole symmetric for \i = 0. 

4-1. Phase diagram in the unbiased case: A = 

The T = mean field phase diagram for a plaquette-centered impurity embedded in 
'intrinsic' (fi = 0) bilayer graphene with A = is shown in figure (3ja). The phase 
diagram shares some qualitative features with that of local moment formation in a 
typical host metal. Namely, there exists a critical ratio of Tq/U before the onset 
of mean field magnetization and a clear Coulomb staircase in the small Tq/U limit. 
Despite these similarities, there are two unusual aspects to this phase diagram. We 
next start by highlighting these novel features and then clarify their physical origin. 

(i) As seen from figure [3ja), there is an extreme skewing of the magnetic regime 
from being centered at (/i — e<j)/Z7 ~ 0.5 for small Tq/U to being centered around 
large positive values of (fi — ed)/U with increasing Tq/U. This strong particle-hole 
asymmetry arises from the fact that the impurity couples asymmetrically to the two 
layers of BLG, leading to a significant real part of the impurity self-energy T,' d {oS). 
The effect of which is to strongly renormalize ej, which causes the observed skewing. 
In order to eliminate this large skewing in later plots, we split the real part of the 
impurity self-energy as 

E» = S' d (0) + (^(u;)-E' d (0)) (17) 

and absorb Sj(0) into the impurity energy, defining a renormalized impurity energy 
£d = (d + E'(0). The resulting renormalized self-energy S^(w) = (T,' d (w) — E^(0)) then 
vanishes at u — 0, and remains small but nonzero away from oj = 0. Plotting the 
impurity model phase diagram in terms of the renormalized impurity energy e^, to a 
large degree but not completely, removes the strong particle-hole asymmetry for fi = 0; 
this can be seen in figure p[b). Of course, strong particle-hole asymmetry continues to 
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Figure 3. Phase diagram of local moment formation on plaquette centered 
impurities in terms of for (a) fi = 0, and in terms of e<2 = e d + (0) for 
(b) fj, = 0, (c) n = -0.05t, and (d) p, = 0.05t d). X = 0.3i in all figures. Grey- 
scale measures the local moment m^. 



exist away from /j, = even after accounting for the impurity energy renormalization, 
as shown in figure [3jc),(d); this can be ascribed to the particle-hole asymmetry of the 
BLG dispersion at nonzero fi. 

(ii) As seen from figure |3ja), there is a dramatic elongation of the magnetic 
region to large values of To/U ~ 10, which one can partially attribute the small 
density of states at /i = 0. However, the phase diagram is also influenced by the 
wavefunctions of the BLG quasiparticles. A close inspection of the phase boundaries 
reveals that they are not symmetric about \x = even after accounting for the self- 
energy correction discussed above. We understand that this residual particle-hole 
symmetry breaking arises from the asymmetric broadening of impurity level caused by 
the disparate effective hybridizations with the different bands. This effect is also seen 
in the phase diagrams for systems by comparing the /i = 0.05t and /i = — 0.05i phase 
diagrams. While one might naively expect that the symmetry between the valence 
and conduction dispersions would lead to symmetric phase diagrams for positive and 
negative chemical potential, subtle di fferences between the two regions again reflect 
the influence of the wavefunctions of the electrons that hybridize with the impurity 
level. 

It is, in fact, extremely instructive to compare the complete impurity phase 
diagram of bilayer graphene with a fictitious system of electrons obtained by setting 
<?iks(r) = exp(ik • r)/-\/Ng in equation |4j where N s is the total number of sites in the 
bilayer. These fictitious electrons are chosen to have the same dispersion as the BLG 
quasiparticles, but their coupling to the impurity does not account for the chirality 
or the band dependence of the quasiparticle wavefunctions. We find that some of the 
unusual features of the BLG impurity phase diagram, discussed above, are eliminated 
upon making this change. 

Most noticeably, the phase diagram of the fictitious fermions is not skewed when 
H = even when plotted in terms of the unrenormalized impurity energy, indicating 
that €d — £d, so that £<j(0) = 0. This stems from a symmetry Ed(— uj) = — Y,* d {iJ) in 
the expression for the self-energy in equation [14] upon assuming a band-independent 
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Figure 4. The phase diagram of local moment formation for fictitious fermions 
with the same dispersion as bilayer graphene for (above) fi = 0, and (below) 
fi = — 0.05i, plotted in terms of e^. (Note, = when fi = 0.) The phase 
diagram for fj, = 0.05t is related to that of fj, = —0.05 by a reflection about 
fi-e d /U = 0.5. 

Vks- Moreover, it also follows that the phase diagram for the fictitious fermions is 
exactly particle-hole symmetric, in contrast to the case of BLG. 

Similar arguments also exactly relate the non-chiral phase diagrams of systems 
with corresponding chemical potential \i and —fi by noting that the self-energy at 
finite chemical potential can be obtained by £<z(^ + A*) of the self-energy at /i = 0. 
Consequently, the phase diagram of the —\i system is obtained by reflecting the phase 
diagram of the fj, system. This is again in contrast to the phase diagram of BLG 
where there is no such relation between systems with positive and negative \i. In BLG, 
particle-hole excitations in a system with positive chemical potential favour different 
bands than those of a system with negative chemical potential. Since each band has 
a unique effective coupling to the impurity in BLG, the hybridization of the impurity 
states will depend on the sign of the chemical potential and so the phase diagrams will 
be different. Finally, the other major distinction between the finite chemical potential 
phase diagrams of the two systems is that, once again, the BLG phase diagram is more 
strongly skewed, even when plotted in terms of E4. This confirms that the band- and 
momentum-dependence of the hybridization to the BLG quasiparticles is responsible 
for sizeable shift in the impurity energy via a sizeable real self-energy. 

4-2. Phase diagram in the biased case: A ^ 

We now turn our attention towards a BLG system in a dual-gate configuration. This 
setup allows one to continuously tune the layer bias and the average chemical potential 
independently by applying an external electric field perpendicular to the layers. In the 
presence of a symmetric interlayer bias, the chemical potential remains fixed while a 
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Figure 5. Phase diagram of local moment formation in plaquette centered 
impurities with A = O.Ot, 0.2t, —0.2t as a function of = + (uj = 0, A = 0). 
In both figures \ = 0.34 and fi = 0. 

band gap opens in the bulk electronic spectrum of BLG. In the context of local moment 
formation, this modification to the density of states is expected to substantially change 
the extent to which an impurity state hybridizes with the BLG electrons. In addition 
to this, the impurity energy levels also shift up or down depending on the potential 
of the layer in which it resides. This remarkable ability to alter the energy of an 
impurity level with respect to the chemical potential through the application of an 
external electric field is unique to multilayer systems, and has no analog in monolayer 
graphene. 

In the first part of this section, we explore how biasing the layers affects local 
moment formation by reconstucting phase diagrams similar to those above, but for 
gated systems with different layer bias and fixed fj, = 0. Doing so allows us to identify 
regions of impurity parameters where local moment formation can be turned on and/or 
off by the electric field. In the subsequent part of this section, we consider the ability to 
tune both the chemical potential and bias by constructing alternative phase diagrams 
where /i and U are varied and it is the bare impurity energy which is fixed. This is 
again done for a selection of values for the bias. 

4.2.1. Impurity energy variation Figure [5] is the phase diagram of the impurity model 
for experimentally accessible values of A plotted in terms of the redefined impurity 
energy ej, introduced above (x = 0.3t and /i = 0). The bias has two effects on the 
impurity model: (i) it opens a band gap ~ A in the BLG dispersion, and (ii) it shifts 
the impurity energy by A/2. Let us discuss, in turn, the impact of these two effects 
on the phase diagram. 

(i) First consider restricting the effect of turning on a bias to opening a gap in 
the BLG spectrum so that the impurity energy level remains unaltered. Then, the 
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dominant effect of a large A is the elongation of the phase boundary to large Tq/U , 
regardless of the parity of A. This occurs because the bias induces a large band 
gap and, when the impurity spectral peak lies in this gap, the coupling between the 
impurity and the extended states becomes negligible because the density of states 
vanishes. (We have to be careful that the renormalized impurity energy, taking self- 
energy corrections into account, should lie in the gap; this renormalization is small 
if the band gap is large compared to IV) Hence, the impurity spectral functions 
become simple delta functions and if we vary for fixed T /U the local moment 
phase boundary resembles that of a simple Coulomb staircase in the atomic limit. 

(ii) The effect of shifting the energy of the impurity level is similar to the effect 
of the real part of the self-energy in the phase diagram; it dramatically skews the local 
moment phase about ed/U = 0.5. The direction of the skewing depends on the parity 
of the bias, as this determines the direction of the impurity energy shift. 

If the impurity energy shift and opening of a band gap are taken together, 
both skewing and elongation of the local moment phase boundary occur. As the 
electric field is increased from zero to large field strengths, the local moment phase 
continuously elongates and 'peels' away from the zero bias boundary. Although slight, 
it is important to note that the phase diagrams with opposite bias parity are not 
symmetric but have slight differences that arise from the breaking of layer symmetry 
by the impurity. One of the key new results is the identification of regions in the 
impurity parameter space where local moments can be turned either on and/or off 
by adjusting the electric field. The region where local moments survive both in the 
presence and absence of the electric field are simply where the phases overlap. 

4-2.2. Chemical potential variation Now we explore the possibility of tuning the 
chemical potential of the system to control local moment formation both in the 
unbiased and biased cases. To do this, we construct phase diagrams for a given 
€d and A, and we now vary /i and U . We do this for A = 0, ±0.2t, for a choice of the 
bare impurity energy such that the noninteracting impurity spectral peak appears in 
the midgap when A = — 0.2t, which we do by choosing + A/2 = — = 0, A). 

In figure [6j the phase diagram is plotted in terms of a redefined impurity energy 
£d = (d + S(cj — 0, A = 0). It is important to emphasize that the location of the 
spectral peaks mostly do not correspond to Ed- The real part of the self-energy has 
significant frequency dependence that shifts the location of the spectral peak, whose 
effect must also be accounted for in order to fully understand the phase diagrams. 

When the system is unbiased (i.e. A = 0) the impurity energy level lies within the 
conduction bands (see reference [T| or reference [23J for details on the band structure) . 
In this case, the phase diagram is qualitatively similar to that of a single site impurity 
(see reference [34 ). When a positive bias is in place, A = Q.2t, a band gap opens 
and the impurity energy shifts deeper into the conduction band. Consequently, the 
phase boundaries for local moment formation are significantly reduced because of the 
enhanced broadening due to the increase in the density of states at higher energy in 
the conduction bands. 

The more interesting case is when A = — 0.2t and the impurity spectral peak 
shifts down into the middle of the gap. Then, if U is small enough so that the doubly 
occupied state also lies at sub-gap energies, both the singly and doubly occupied 
states can no longer hybridize with the BLG states and the impurity spectral function 
reduces to delta functions. Hence, we again recover local moment formation very 
similar to the atomic limit, but now in the large Tq/U limit. However, in this limit 
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Figure 6. Local moment phase diagram for biased bilayer graphene. The 
impurity level at A = — 0.2t bias was chosen so that its spectral peak lies in 
the middle of the gap. 



the upper and lower phase boundaries of the Coulomb-staircase are not separated by 
(/i — £d)/U — 1 because of level repulsion. The doubly occupied state shifts down in 
energy due to T,' d ((jj). 

Thus, this phase diagram is unusual in the sense that it has two regimes resembling 
the atomic limit at large and small Tq/U. Separating these regimes is the part of the 
phase diagram where the doubly occupied state's energy lies beyond the band edge 
and hybridizes with the conduction states. This occurs at about T /U ~ 2.5 when 
U ~ O.Lt, precisely where the unusual 'hump'-like feature is seen in the upper phase 
boundary. The cause of the feature can again be attributed to level repulsion, as it 
becomes very strong for states close to the the gap edge and S'(cj) exhibits a large 
peak. 

5. RKKY interaction between local moments 

In this Section, we explore the RKKY coupling between local moments [37J EH1 EH) 
and study how it can be tuned by varying the band gap and chemical potential using 
a dual-gate configuration. We have seen in the previous section that such variations 
will, in general, modify the local moment. Here, we focus on changes to the RKKY 
coupling induced purely by changes in the bulk band structure and filling. 

We consider two classical local moments that couple to the set of sites {r} and 
{r'}, respectively, 



where Jr is the strength of the exchange coupling of an electron's spin, s r , at 
site r with the magnetic impurity S a . Upon integrating out the itinerant electrons 
and retaining only those terms that are second order in J* , one obtains a reduced 
Hamiltonian for the local moments, 




(18) 



Heff — JRKKYSl ■ S2. 



(19) 
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Figure 7. Crystal structure of bilayer graphene and site labelling convention. 
The primitive lattice vectors are a and b and the armchair and zigzag directions 
are indicated by the arrows. The local moments considered here are plaquette 
centered and reside on the same layer above an A2 atom. 



The coupling Jrkky is given by 

Jrkky=^ £ M^(q)0rW^(i)^ q (i)C +q (Oe iq - (ri - r2) 

qki j nm 

x M ^'7 (S) . <») 

^k+q <?k 

where m/n are band indices, i/j are the combined sublattice and layer label, np is the 
Fermi distribution, and My(q) is a matrix describing the Fourier transform between 

different sites weighted by Jr'J^,. The explicit form of My for the case of interest 
is provided below. 

For monolayer graphene, it has been shown that a perturbative treatment in the 
continuum low-energy theory [TJ] produces approximate results that match closely 
with exact diagonalization [Tl] and lattice Green's functions methods [TBJ, as long 
as an appropriate high-energy cutoff scheme is applied. In the above perturbative 
treatment, the entire band structure is used in the calculation so as to avoid any 
cutoff dependence and the RKKY coupling is accurately reproduced for monolayer 
graphene. We therefore expect this perturbative calculation to also be a reasonable 
approach to study the RKKY coupling in BLG in the dual-gate configuration. 

We analyzed various moment configurations such as single site (AA, BB, AB) and 
plaquette coupled moments both along the zigzag and armchair directions (see figure 
for labelling conventions) . The effects of varying the chemical potential and layer bias 
were seen to be qualitatively similar for each case, so we have chosen to present only 
the results for plaquette centered moments that lie along the zigzag direction. The 
impurity atom is taken to lie above an A2 in the center of a hexagonal plaquette in 
layer 1. For simplicity, we assume the coupling to each of the seven sites is equal 
so that Jr 1 ' 1 = Jr 2 ^ = J, although we have checked that the results are qualitatively 
unaffected if there is an unequal coupling to the site below the impurity on the other 
layer (A2). In this case, the components of M$j are 

M AlAl = J 2 (4 + 2 (cos(g a ) + cos((? b ) + cos(g a - q b )) J 
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M AlBl = J 2 (2 + 2 (cos(q a ) + co S (q b ) + e *(«.-*) + e -«39.-»)) 
M AlA2 = J 2 (l + e lqb + e^-^A 
M BlA . 2 = J 2 (e lq " + e ^- qb) + e l{q "- 2qb) ^ 

M A2 a 2 =J 2 , M,, .u;,. (21) 

where q a — q • a, q b = q • b and a = x and b = x/2 + v3y/2. 

To demonstrate the ability to tune the RKKY interaction using the dual-gate 
configuration, the Jrkky coupling, normalized to its value at fi = and A = 0, is 
plotted in figure [H] as a function of the interlayer bias A for two moments separated 
by 10 lattice spacings. This is done for \x = and /1 = 0.054. For experimental 
considerations, one must keep in mind that the RKKY coupling is quite small in 
BLG. As an example, a bare exchange term equal to JW = — 0.24 produces 
an effective coupling Jrkky — 1-3 x 10 _4 4 (~ 4.4 K) at 4 lattice spacings, and 
just Jrkky = 7.5 x 10~ 6 4 (~ 0.3 K) at 10 lattice spacings. However, similar to 
monolayer graphene, electron interactions are expected to make the coupling strength 
more long ranged [15] . At shorter distances, the RKKY interaction is enhanced, but 
the tunability is reduced. 

Before describing the tunable features of the RKKY coupling, it is important to 
first understand that the wavefunctions of a given band are sensitive to the parity 
of the bias between the layers, even though the dispersion is not. Their dependancy 
on the parity can significantly influence how Jrkky changes with bias, as explained 
below. When a positive bias is present, states in the upper two bands are more heavily 
weighted to layer 1 sites, while states in the lower band are more heavily weighted to 
the layer 2 sites. This weighting is reversed when the parity of the bias is negative. 
In contrast, when there is no bias the weighting of the wavefunction is the same for 
each layer. 

With this background, it is possible to explain the symmetry/asymmetry between 
the two curves. When fj, = 0, the chemical potential lies between the valence and 
conduction bands, and so particle-hole excitations can only occur between them. This 
corresponds to one of the states being localized to layer 1 and the other localized to 
layer 2. It follows that the coupling strength Jrkky is parity invariant and so it is 
symmetric for positive and negative biases. 

In contrast, when \i ^ 0, the coupling is sensitive to the parity of the bias. 
If = 0.054, A > and \i is less than the band gap, the chemical potential lies 
in the third band where the states tend to localize to layer 1, the layer in which 
the moments reside. The finite chemical potential causes some of the particle-hole 
excitations between the lower and upper two bands to be suppressed by Pauli-blocking, 
and also introduces low-energy excitations between the two upper bands where the 
wavefunctions are weighted to layer 1. If however, fj, = 0.054, A < and /i is less than 
the band gap, the chemical potential lies in the third band, but now these states tend 
to localize to layer 2. Although the energetics of the scattering processes remain the 
same, the matrix elements do not. The excitations between the upper two bands now 
have matrix elements whose weighting on layer 1 is much less. Thus, the coupling 
is dependent on the relative parity o f /i to A. Hence, if wc consider a system with 
fi = —0.054, the Jrkky curve will be reflected about A = 0. 

In addition to the effects described above, the density of states about the chemical 
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Figure 8. Normalized RKKY coupling strength between to classical plaquette 
centered moments at a distance of 10 lattice spacings along the zigzag direction. 
Both moments are located on the same layer and the chemical potential is chosen 
to be fx = and fj, = 0.05i at a temperature T = 0.0024 (~ 70 K). 



potential tends to increase for small biases, as the band edge flattens and is pushed 
closer to the chemical potential. At large bias strengths, the dispersion close to the 
band-edge resembles that of a 'mexican-hat', leading to further complexity in the 
density of states. Furthermore, at finite chemical potential, the Fermi points extend 
out to form a Fermi-surface symmetric about the if-points. The combination of all 
these effects lead to the non-trivial changes seen the RKKY coupling in figure [Hj 

Interestingly, at a distance of 10 lattice spaces, the coupling remains 
antiferromagnetic when /i = 0.05i and A > and tends to increase with bias 
strength. However, when A becomes increasingly negative, the antiferromagnetic 
coupling strength is reduced to zero then switches to a ferromagnetic coupling about 
A ~ —0.12. For the case of /i = 0.05i considered here, the ability to fully turn 
off the coupling and switch the sign of Jrkky with an applied electric field sets 
in when the moments are separated by at least 9 lattice spacings and persists to 
about 20 lattice spacings. This window of tunability may be augmented by carefully 
adjusting \i. Regardless of the sign of the bias, once the band gap exceeds the chemical 
potential, the /i = 0.05i curve begins to merge with the /i = curve, as expected at 
low temperature. 

This ability to dynamically tune the strength of the RKKY interaction, as well as 
its sign from being antiferromagnetic to being ferromagnetic, for two fixed moments 
using just an electric field (rather than doping) is perhaps the most interesting feature 
of this system. 

6. Discussion 

In this paper, we have discussed the physics of local moments on bilayer graphene. 
We have, in our discussion, ignored the effects of electron-electron interactions among 
the BLG electrons. These are known to be important for quadratic band touching 
[40l 1411 l42l l43l [44] , leading to symmetry breaking and many-body gaps for zero doping 
and zero electric field. However, so long as there is nonzero doping or the presence of 
an electric field that gaps out the low-energy BLG states, the perturbative effects of 
electron-electron interactions are benign and will not lead to qualitatively new many- 
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body effects. Using a clean substrate may also be a viable route to mitigating the 
effects of rippling and disorder. Moreover, the substrate may partially screen the 
BLG electron-electron interactions and reduce the many-body gap recently observed 
in suspended BLG [45] . The competition between impurity physics and many-body 
interactions in BLG deserves a careful separate investigation. The ability to turn 
on/off local moments, and the ability to tune the sign and magnitude of the RKKY 
coupling between local moments using electric fields perpendicular to the bilayer, 
which we have studied, constitutes physics beyond what has been discussed for 
monolayer graphene. The experimental realization of such tunable local moments 
in BLG is a compelling prospect. It would be interesting to study such a system using 
scanning tunnelling spectroscopy and to probe the quantum dynamics of interacting 
local moments in experiments. We expect that thermal fluctuations will only slightly 
alter the phase boundaries as long as the temperature is below the Hubbard gap. 
Quantum fluctuations are expected to lead to Kondo screening or valence fluctuations 
in points of the phase diagram - this is an interesting direction for future research. 
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